#pragma GCC optimize(3)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <pthread.h>
#include <memory.h>

unsigned long long thread_count = 15;

// M : number of combinations to store
const int M = 1000;
// dim : dimension of metric space
int dim;
// n : number of points
int n;
// k : number of pivots
int k;

long long total;

double *coord;
double *maxDistanceSum;
int *maxDisSumPivots;
double *minDistanceSum;
int *minDisSumPivots;

pthread_t *thread_handles;
pthread_spinlock_t sort_spinlock;

long long fact(int n)
{
    long long factorial = 1;
    for (; n >= 1; n--)
    {
        factorial *= n;
    }
    return factorial;
}

long long com(int n, int m)
{
    long long result = 1L;
    int temp_num = n;
    for (int i = 0; i < m; ++i)
    {
        result *= temp_num;
        temp_num--;
    }
    result /= fact(m);

    return result;
}

// Calculate sum of distance while combining different pivots. Complexity : O( n^2 )
double sumDistance(int *pivots, double **max_table, double *rebuiltCoord, int use_max_table)
{
    double chebyshevSum = 0;
    if (use_max_table)
        for (int i = 0; i < n; i++)
            for (int j = i; j < n; j++)
            {
                double chebyshev = max_table[i][j];
                double dis = fabs(rebuiltCoord[i * k + k - 1] - rebuiltCoord[j * k + k - 1]);

                chebyshev = dis > chebyshev ? dis : chebyshev;
                chebyshevSum += 2 * chebyshev;
            }

    else
        for (int i = 0; i < n; i++)
            for (int j = i; j < n; j++)
            {
                double chebyshev = 0;
                int ki;
                for (ki = 0; ki < k; ki++)
                {
                    double dis = fabs(rebuiltCoord[i * k + ki] - rebuiltCoord[j * k + ki]);
                    chebyshev = dis > chebyshev ? dis : chebyshev;
                }
                chebyshevSum += 2 * chebyshev;
            }

    return chebyshevSum;
}

void calc(int *pivots, double **max_table, double *rebuiltCoord, int use_max_table)
{
    // Calculate sum of distance while combining different pivots.
    double distanceSum = sumDistance(pivots, max_table, rebuiltCoord, use_max_table);

    pthread_spin_lock(&sort_spinlock);
    // put data at the end of array
    maxDistanceSum[M] = distanceSum;
    minDistanceSum[M] = distanceSum;
    int kj;
    for (kj = 0; kj < k; kj++)
    {
        maxDisSumPivots[M * k + kj] = pivots[kj];
    }
    for (kj = 0; kj < k; kj++)
    {
        minDisSumPivots[M * k + kj] = pivots[kj];
    }
    // sort
    int a;
    for (a = M; a > 0; a--)
    {
        if (maxDistanceSum[a] > maxDistanceSum[a - 1])
        {
            double temp = maxDistanceSum[a];
            maxDistanceSum[a] = maxDistanceSum[a - 1];
            maxDistanceSum[a - 1] = temp;
            int kj;
            for (kj = 0; kj < k; kj++)
            {
                int temp = maxDisSumPivots[a * k + kj];
                maxDisSumPivots[a * k + kj] = maxDisSumPivots[(a - 1) * k + kj];
                maxDisSumPivots[(a - 1) * k + kj] = temp;
            }
        }
    }
    for (a = M; a > 0; a--)
    {
        if (minDistanceSum[a] < minDistanceSum[a - 1])
        {
            double temp = minDistanceSum[a];
            minDistanceSum[a] = minDistanceSum[a - 1];
            minDistanceSum[a - 1] = temp;
            int kj;
            for (kj = 0; kj < k; kj++)
            {
                int temp = minDisSumPivots[a * k + kj];
                minDisSumPivots[a * k + kj] = minDisSumPivots[(a - 1) * k + kj];
                minDisSumPivots[(a - 1) * k + kj] = temp;
            }
        }
    }

    pthread_spin_unlock(&sort_spinlock);

    // if (pivots[k - 1] == n - 1)
    // {
    //     for (int i = 0; i < k; ++i)
    //         printf("%d ", pivots[i]);
    //     printf("finished\n");
    // }

    return;
}

void *CALCULATOR(void *rank)
{
    printf("CALCULATOR线程%d开始运行\n", rank);

    long long my_rank = (long long)rank;

    long long floor, ceil; // 范围 [floor, ceil)

    long long index = 0;

    long long per = total / thread_count; // 每个线程的工作量

    int *temp = (int *)malloc(sizeof(int) * k);

    int *max_temp = (int *)malloc(sizeof(int) * (k - 1));

    double **max_table = (double **)malloc(n * sizeof(double *));
    for (int i = 0; i < n; i++)
        max_table[i] = (double *)malloc(n * sizeof(double));
    for (int i = 0; i < n; i++)
        max_table[i][i] = 0;

    double *rebuiltCoord = (double *)malloc(sizeof(double) * n * k);
    for (int i = 0; i < n * k; i++)
        rebuiltCoord[i] = 0;

    floor = my_rank * per;
    ceil = (my_rank == thread_count - 1) ? total : (my_rank + 1) * per;

    for (int i = 0; i < k; i++)
        temp[i] = i; // 初始化为0, 1, 2, 3, ...

    for (int i = 0; i < k - 1; i++)
        max_temp[i] = -1;

    while (index != ceil)
    {
        if (index >= floor)
        {
            // Rebuild coordinates. New coordinate of one point is its distance to each pivot.
            for (int i = 0; i < n; i++)
            {
                for (int ki = 0; ki < k; ki++)
                {
                    double distance = 0;
                    int pivoti = temp[ki];
                    for (int j = 0; j < dim; j++)
                        distance += pow(coord[pivoti * dim + j] - coord[i * dim + j], 2);
                    rebuiltCoord[i * k + ki] = sqrt(distance);
                }
            }


            if (k > 2)
            {
                int changed = 0;
                for (int i = 0; i < k - 1; i++)
                {
                    // if (temp[i] != max_temp[i] && index < (ceil - floor) / boundary)
                    if (temp[i] != max_temp[i])
                        changed = 1;
                }
                if (changed)
                {
                    for (int i = 0; i < n; i++)
                        for (int j = i; j < n; j++)
                        {
                            double max_num = 0;
                            for (int ki = 0; ki < k - 1; ki++)
                            {
                                // double dis = fabs(distTable[i][temp[ki]] - distTable[j][temp[ki]]);
                                double dis = fabs(rebuiltCoord[i * k + ki] - rebuiltCoord[j * k + ki]);
                                max_num = dis > max_num ? dis : max_num;
                            }
                            max_table[i][j] = max_num;
                            max_table[j][i] = max_num;
                        }

                    for (int i = 0; i < k - 1; i++)
                        max_temp[i] = temp[i];
                }
                calc(temp, max_table, rebuiltCoord, 1);
            }
            else
                calc(temp, max_table, rebuiltCoord, 0);
        }

        index++;
        temp[k - 1]++;
        for (int i = k - 1; i > 0; --i)
        {
            if (temp[i] == n - k + i + 1)
            {
                temp[i - 1]++;
                for (int j = i; j < k; j++)
                    temp[j] = temp[j - 1] + 1;
            }
        }
    }

    printf("CALCULATOR线程%d结束运行\n", rank);
    free(temp);
    free(max_temp);
    free(max_table);
    free(rebuiltCoord);
    return NULL;
}

int main(int argc, char *argv[])
{
    // filename : input file namespace
    char *filename = (char *)"uniformvector-2dim-5h.txt";

    if (argc == 2)
    {
        thread_count = strtol(argv[1], NULL, 10);
    }
    else if (argc == 3)
    {
        thread_count = strtol(argv[1], NULL, 10);
        filename = argv[2];
    }
    else if (argc != 1)
    {
        printf("Usage: ./pivot <thread_count> <filename>\n");
        return -1;
    }

    thread_handles = (pthread_t *)malloc(thread_count * sizeof(pthread_t));

    // Read parameter
    FILE *file = fopen(filename, "r");
    if (file == NULL)
    {
        printf("%s file not found.\n", filename);
        return -1;
    }
    fscanf(file, "%d", &dim);
    fscanf(file, "%d", &n);
    fscanf(file, "%d", &k);
    printf("dim = %d, n = %d, k = %d\n", dim, n, k);

    // Start timing
    struct timeval start;

    // Read Data
    coord = (double *)malloc(sizeof(double) * dim * n);
    int i;
    for (i = 0; i < n; i++)
    {
        int j;
        for (j = 0; j < dim; j++)
        {
            fscanf(file, "%lf", &coord[i * dim + j]);
        }
    }
    fclose(file);
    gettimeofday(&start, NULL);

    // maxDistanceSum : the largest M distance sum
    maxDistanceSum = (double *)malloc(sizeof(double) * (M + 1));
    for (i = 0; i < M; i++)
    {
        maxDistanceSum[i] = 0;
    }
    // maxDisSumPivots : the top M pivots combinations
    maxDisSumPivots = (int *)malloc(sizeof(int) * k * (M + 1));
    for (i = 0; i < M; i++)
    {
        int ki;
        for (ki = 0; ki < k; ki++)
        {
            maxDisSumPivots[i * k + ki] = 0;
        }
    }
    // minDistanceSum : the smallest M distance sum
    minDistanceSum = (double *)malloc(sizeof(double) * (M + 1));
    for (i = 0; i < M; i++)
    {
        minDistanceSum[i] = __DBL_MAX__;
    }
    // minDisSumPivots : the bottom M pivots combinations
    minDisSumPivots = (int *)malloc(sizeof(int) * k * (M + 1));
    for (i = 0; i < M; i++)
    {
        int ki;
        for (ki = 0; ki < k; ki++)
        {
            minDisSumPivots[i * k + ki] = 0;
        }
    }

    total = com(n, k);

    pthread_spin_init(&sort_spinlock, PTHREAD_PROCESS_PRIVATE);
    long long thread;

    for (thread = 0; thread < thread_count; ++thread)
        pthread_create(&thread_handles[thread], NULL, CALCULATOR, (void *)thread);
    for (thread = 0; thread < thread_count; ++thread)
        pthread_join(thread_handles[thread], NULL);

    pthread_spin_destroy(&sort_spinlock);

    // End timing
    struct timeval end;
    gettimeofday(&end, NULL);
    printf("Using time : %f ms\n", (end.tv_sec - start.tv_sec) * 1000.0 + (end.tv_usec - start.tv_usec) / 1000.0);

    // Store the result
    FILE *out = fopen("result-pthread-big.txt", "w");
    for (i = 0; i < M; i++)
    {
        int ki;
        for (ki = 0; ki < k - 1; ki++)
        {
            fprintf(out, "%d ", maxDisSumPivots[i * k + ki]);
        }
        fprintf(out, "%d\n", maxDisSumPivots[i * k + k - 1]);
    }
    for (i = 0; i < M; i++)
    {
        int ki;
        for (ki = 0; ki < k - 1; ki++)
        {
            fprintf(out, "%d ", minDisSumPivots[i * k + ki]);
        }
        fprintf(out, "%d\n", minDisSumPivots[i * k + k - 1]);
    }
    fclose(out);

    // Log
    int ki;
    printf("max : ");
    for (ki = 0; ki < k; ki++)
    {
        printf("%d ", maxDisSumPivots[ki]);
    }
    printf("%lf\n", maxDistanceSum[0]);
    printf("min : ");
    for (ki = 0; ki < k; ki++)
    {
        printf("%d ", minDisSumPivots[ki]);
    }
    printf("%lf\n", minDistanceSum[0]);

    free(thread_handles);
    free(coord);
    free(maxDistanceSum);
    free(maxDisSumPivots);
    free(minDistanceSum);
    free(minDisSumPivots);

    system("pause");

    return 0;
}
